The present invention relates generally to nuclear magnetic resonance imaging methods and systems, and in particular, relates to segmentation of a human internal organ, or a portion of an internal organ, for example a left ventricular epicardium.
When a substance such as human tissue is subjected to a uniform magnetic field (polarizing field B0), the individual magnetic moments of the spins in the tissue attempt to align with this polarizing field, but precess about it in random order at their characteristic Larmor frequency. If the substance, or tissue, is subjected to a magnetic field (excitation field B1) which is in the x-y plane and which is near the Larmor frequency, the net aligned moment, Mz, may be rotated, or xe2x80x9ctippedxe2x80x9d, into the x-y plane to produce a net transverse magnetic moment Mt. A signal is emitted by the excited spins after the excitation signal B1 is terminated, this signal may be received and processed to form an image.
When utilizing these signals to produce images, magnetic field gradients (Gx Gy and Gz) are employed. Typically, the region to be imaged is scanned by a sequence of measurement cycles in which these gradients vary according to the particular localization method being used. The resulting set of received NMR signals are digitized and processed to reconstruct the image using one of many well-known reconstruction techniques.
Most NMR scans currently used to produce medical images require many minutes to acquire the necessary data. The reduction of this scan time is an important consideration since reduced scan time increases patient throughput, improves patient comfort, and improves image quality by reducing motion artifacts. There is a class of pulse sequences which have a very short repetition time (TR) and result in complete scans which can be conducted in seconds rather than minutes. When applied to cardiac imaging, for example, a complete scan from which a series of images showing the heart at different phases of its cycle can be acquired in a single breath-hold.
The prognosis of patients with a wide variety of cardiac diseases (including coronary artery disease, valvular heart disease, congestive heart failure and cardiac arrhythmias) has been closely linked to the performance of the heart as indicated by measurements such as wall thickening, wall motion, and myocardial mass. Accurate quantitative measures of regional contractile function could therefore have significant prognostic and therapeutic importance. For example, many patients with severe coronary artery disease may have normal regional and global left ventricular function at rest but have abnormalities induced by stress. In clinical practice, patients with coronary artery disease can be detected by stress echocardiography based on new functional deficits during stress. However, interobserver variability of this type of qualitative measure is an inherent limitation that could be improved with quantitative measures. Thus, there is a need for high quality quantitative measures of regional cardiac function.
Segmentation of the left ventricle in MR images is therefore a fundamental step in analyzing the performance of the heart. MR image data of the epicardial boundary is currently acquired by applying a specific sequence of RF pulses to yield a NMR signal that provides information pertaining to the tissue under test. A particular pulse sequence can therefore be applied to obtain a mask of pixels in the intensity range of, for example, a cross-section of the left ventricle tissue. Current processes are available for segmenting the epicardium, but they lack robustness and are difficult to use.
Segmentation methods that are currently available include snake-based techniques such as that described by A. Yezzi, et al. xe2x80x9cA Geometric Snake Model for Segmentation of Medical Imagery,xe2x80x9d IEEE Transaction on Medical Imaging, 16, 199-209 (April, 1997). Snakes, also known as active contours, have been used in an attempt to segment features of the left ventricle. Snakes are described by a parameterized curve whose evolution is determined by the minimization of an energy field. The equation of the energy field, as defined by J. C. Gardner et al. xe2x80x9cA Semi-Automated Computerized System for Fracture Assessment of Spinal X-Ray Films,xe2x80x9d Proceedings of the International Society for Optical Engineering, 2710, 996-1008 (1996), is:                               E          ⁡                      [                                          x                →                            ⁢                              xe2x80x83                            ⁢                              (                s                )                                      ]                          ≡                  k          ⁢                      xe2x80x83                    ⁢                                    ∫              0              1                        ⁢                          xe2x80x83                        ⁢                          ⅆ                              s                ⁡                                  [                                                                                    1                        2                                            ⁢                                              xe2x80x83                                            ⁢                      α                      ⁢                                              xe2x80x83                                            ⁢                                                                        (                                                                                    ⅆ                                                              x                                →                                                                                                                    ⅆ                              s                                                                                )                                                2                                                              +                                                                  1                        2                                            ⁢                                              xe2x80x83                                            ⁢                      β                      ⁢                                              xe2x80x83                                            ⁢                                                                        (                                                                                    ⅆ                                                                                                                                   2                                                                ⁢                                                                  x                                  →                                                                                                                                                    ⅆ                                                              s                                2                                                                                                              )                                                2                                                              -                                          γ                      ⁢                                              xe2x80x83                                            ⁢                      H                      ⁢                                              xe2x80x83                                            ⁢                                              (                                                                              x                            →                                                    ⁢                                                      xe2x80x83                                                    ⁢                                                      (                            s                            )                                                                          )                                                                              ]                                                                                        (        1        )            
where s is the parameterization variable, {overscore (x)} is the parameterized curve, xcexa is the normalization constant, xcex1 is the H({overscore (x)})=|{overscore (∇)}/({overscore (x)})| tension of the snake, xcex2 is the rigidity of the snake, xcex3 controls the attraction to image features, and I is the pixel intensity of the image. H(x) refers to a function which defines the features that attract the snake algorithm to the boundary and, typically, is chosen to be the magnitude of the gradient of the image intensity.
Because the magnitude of the gradient is used to attract the algorithm to the boundary of the left ventricle, the snake does not work well where the boundary is defined by edges that are weak in intensity. In order for the snake algorithm to attach to a boundary, a user must intervene and supply a boundary condition to define the proximity of the boundary for the snake. This is undesirable because the user may need to interact with the segmentation algorithm while the images are being processed. Snake based techniques can be used, as described by Yezzi, to produce a geometric snake model having a stopping term and a constant inflation term added to the evolution equation. The resulting evolution equation of the Yezzi active contour model is:                                           ∂            Ψ                                ∂            t                          =                              φ            ⁢                          xe2x80x83                        ⁢                          "LeftDoubleBracketingBar"                              ∇                Ψ                            "RightDoubleBracketingBar"                        ⁢                          xe2x80x83                        ⁢                          (                              κ                +                v                            )                                +                                    ∇              φ                        *                          ∇              Ψ                                                          (        2        )            
here v is a constant inflation force,   κ  ≡      div    ⁢          xe2x80x83        ⁢          (                        ∇          ψ                          "LeftDoubleBracketingBar"                      ∇            ψ                    "RightDoubleBracketingBar"                    )      
is the curvature of the level sets of "psgr"(x, y, t), xcfx86 is a function dependent on the type of image and is a stopping term for the curve evolution. Snake based techniques are additionally unfavorable because they rely primarily on edge information only, and therefore are subject to greater error and generally lack robustness, particularly in a clinical setting. S. Ranganath attempted unsuccessfully to segment an epicardium using a snake, as described in xe2x80x9cContour Extraction from Cardiac MRI Studies Using Snakes,xe2x80x9d IEEE ransactions on Medical Imaging, 14(2), 328-338 (June, 1995).
Another such method currently used in conjunction with attempted detection of epicardial boundaries is a shape-based technique known as the MR Analytical Software System (MASS), introduced by R. J. van der Geest et al. xe2x80x9cComparison Between Manual and Semiautomated Analysis of Left Ventricular Volume Parameters from Short-Axis MR Images,xe2x80x9d Journal of Computer Assisted Tomogrophy,xe2x80x9d 21(5), 756-675 (1997), which uses shape as the central principal for the detection of the epicardial and endocardial contours. The MASS algorithm operated by first using a Hough transform, well known in the art, to determine the initial search location for the endocardial and epicardial boundaries. The Hough transform produces a map with high values near the center of approximately circular objects in the original image. A size constraint is then used to narrow a search for circular areas in the image corresponding to the first cardiac phase. After the search determines which circular areas constitute the boundary areas, a line is fit through the Hough images to estimate the center of the left ventricle. The line provides an estimate of the longitudinal axis of the heart.
The MASS algorithm then transforms each image in the study to a polar image and computes a polar edge image. Using a circle estimation from the original image, the intensity of edges in the radial direction, an estimate for myocardial wall thickness, and a maximum likelihood estimate of the endocardial and epicardial radii are calculated. If a satisfactory estimate is not found for the epicardial radius, one is created afterward through linear interpolation between adjacent radii. Once the epicardial boundary has been determined, MASS uses an intensity thresholding technique to find the endocardial boundary. However, because shape-based techniques primarily rely on the shape of the image to produce the outer edge pattem, these methods, like the snake, are subject to error and generally lack robustness.
What is therefore needed is a method and apparatus for segmenting an epicardium in an image that relies on several information sources to produce an image of the left ventricular epicardial boundary that is clinically robust and that operates with greater accuracy than conventional techniques and that requires only minimal user interaction.
The present invention relates to a system and method for segmenting a human organ, and in particular, a left ventricular epicardium using a method that relies on image shape, gradients, and intensity, and requires only minimal user input to provide a clinically robust mask image of the epicardial boundary of the left ventricle of the heart.
In accordance with a first aspect of the invention, a method for segmenting an image acquired with a medical imaging system to identify the boundary of an organ includes acquiring image data of the organ, and subsequently reconstructing an image of the organ using the acquired image data. Next, an intensity map is produced by segmenting the reconstructed image to produce pixels lying within the boundary of the organ. Next, an edge map is created by detecting the edges of the reconstructed image, and information from the edge map is used to refine the intensity map so as to include variations in intensity corresponding to the outer boundary of the organ. Once the intensity map is refined, a center of the organ is located and the intensity map is searched outwardly from the center to locate the variations in intensity. An image of the outer boundary corresponding to the variations is then generated.